##################################NOTE###################################
Use 120423_Fmask_alpha.R instead of this script.
################################## END ##################################

##################################NOTE###################################
Use 120423_Fmask_alpha.R instead of this script.
################################## END ##################################

##################################NOTE###################################
Use 120423_Fmask_alpha.R instead of this script.
################################## END ##################################

##################################NOTE###################################
Use 120423_Fmask_alpha.R instead of this script.
################################## END ##################################

#TEST#
#' Detect cloud, cloud shadow and snow using the Fmask method
#'
#' The Fmask (Function of mask) is an automated cloud, cloud shadow and snow masking detect ion method for Landsat TM/ETM+ images developed by Zhu and Woodcock (2012). A Matlab code for the Fmask is maintained by Zhu (http://code.google.com/p/fmask/). This is translation of the Matlab code to R. 
#'
#' @name FmaskforR-package
#' @aliases FmaskforR
#' @docType package
#' @title Automated cloud, cloud shadow and snow masking detection method for Landsat TM/ETM+ images.
#' @author Joseph Henry \email{joseph.henry@@sydney.edu.au}, Willem Vervoort \email{willem.vervoort@@sydney.edu.au}, Denis Duro and Mark Adams
#' @references Zhu, Z. and Woodcock, C.E. 2012. Object-based cloud and cloud shadow detection in Landsat imagery. Remote Sensing of Environment 118, 83-94. ' \url{http://code.google.com/p/fmask/}
#' @keywords package, Fmask, Landsat, Cloud, Snow, Shadow


### CONDITION OF USE ###
# Joseph reserves the right to use this code or its derivatives in the creation of an R package.

### ALPHA VERSION ###

## Packages (dependicies of the package)
if (!require(stringr)){
        install.packages('stringr')
		require(stringr)
} 
if (!require(raster)){
        install.packages('raster')
		require(raster)
}
if (!require(RSAGA)){
        install.packages('RSAGA')
		require(RSAGA)
}  
if (!require(rgdal)){
        install.packages('rgdal')
		require(rgdal)
}  
# if (!require(ff)){
        # install.packages('ff')
		# require(ff)
# }  
# if (!require(ffbase)){
        # install.packages('ffbase')
		# require(ffbase)
# }  
if (!require(igraph0)){
        install.packages('igraph0')
		require(igraph0)
}  
# if (!require(bigmemory)){
        # install.packages('bigmemory')
		# require(bigmemory)
# }  

 # require (stringr)
 # require (raster)
 # require(RPyGeo)
	 

# ## ArcGIS setup
# fill_func <- function (in_surface_raster, out_surface_raster, z_limit) {
	# rpygeo.geoprocessor (ARCGIS.FILLSA,
	# c (in_surface_raster, out_surface_raster, z_limit),
	# env = env.cur,
	 # working.directory = reproj_d)
	# }
# ## Configure python
  # env.cur <- rpygeo.build.env(
  # workspace = "C:/",# ARC related
  # python.path = "C:/Python25",python.command = "python.exe" #python related
  # ,overwriteoutput = 1 #if 1, overwrite
  # )
	# # Working directory
   # reproj_d <- paste("C:/",sep="")


# shad_cell_func <- function (input.a, input.b){
# #' Function used in Fmask.
	# shadow_cells <- cellFromRowCol (input.a, input.b[,1], input.b[,2])
	# return (shadow_cells)
#}

#FmaskForR <- function (INPUT_BANDS, INPUT_BANDS_DN, run.prob.over.water, L7.SUN.AZIMUTH.TAN.RADIANS, L7.SUN.AZIMUTH, L7.SUN.ELEVATION.TAN.RADIANS,OUTPUT) {
#' Produce Landsat reflectance and temperature rasters.
#' @param INPUT_BANDS file paths of reflectance and temperature data (bands 1 to 7) (character)
#' @param INPUT_BANDS_DN file paths of digital number data (bands 1 to 3) (character)
#' @param run.prob.over.water choice,  "yes" or  "no", to run Fmask over water. JH didn't use this because no large water bodies in catchment.
#' @param OUTPUT file path of output file (character).
#' @param L7.SUN.AZIMUTH.TAN.RADIANS sun azimuth in radians derived from metadata file that comes with tile (numeric)
#' @param L7.SUN.AZIMUTH sun azimuth in degrees derived from metadata file that comes with tile (numeric)
#' @param L7.SUN.ELEVATION.TAN.RADIANS  sun elevation derived from metadata file that comes with tile (numeric)
 ## Load
#################Constants set for testing 
	#DATA_TR_OUT <- list.files(full.names = TRUE, "C:/Users/Uni/Desktop/NANO/USYD/desktop/Fmask_test/Output/TR/LE70900852001003ASA00") # directory with landsat files for one pass only
	INPUT_BANDS <- DATA_TR_OUT
	INPUT_BANDS_DN <- DATA_PATH_RAW_EXCL
	OUTPUT <- "C:/Users/Uni/Desktop/fmask.tif"
#################
	
	L7.FILLVALUE <- 255
	# Radiances
		Band_1_rad <- raster (INPUT_BANDS [1]) # Band 1
		gc()
		Band_1_rad_vec <- getValues(Band_1_rad) # Vectors
		gc()
		B1.OUT.PIX <- which (Band_1_rad_vec == Band_1_rad@data@min)
		rm (Band_1_rad)
		rm (Band_1_rad_vec)

		Band_2_rad <- raster (INPUT_BANDS [2]) # Band 2
		gc()
		Band_2_rad_vec <- getValues(Band_2_rad) # Vectors		
		gc()
		B2.OUT.PIX <- unique (c (B1.OUT.PIX, which (Band_2_rad_vec == Band_2_rad@data@min)))
		rm (Band_2_rad)
		rm (B1.OUT.PIX )
		rm (Band_2_rad_vec)
		gc()

		Band_3_rad <- raster (INPUT_BANDS [3]) # Band 3
		gc()
		Band_3_rad_vec <- getValues(Band_3_rad) # Vectors		
		gc()
		B3.OUT.PIX <- unique (c (B2.OUT.PIX, which (Band_3_rad_vec == Band_3_rad@data@min)))
		rm (Band_3_rad)
		rm (B2.OUT.PIX )
		rm (Band_3_rad_vec)
		gc()

		Band_4_rad <- raster (INPUT_BANDS [4]) # Band 4
		gc()
		Band_4_rad_vec <- getValues(Band_4_rad) # Vectors	
		gc()
		B4.OUT.PIX <- unique (c (B3.OUT.PIX, which (Band_4_rad_vec == Band_4_rad@data@min)))
		rm (Band_4_rad)
		rm (B3.OUT.PIX )
		rm (Band_4_rad_vec)
		gc()

		Band_5_rad <- raster (INPUT_BANDS [5]) # Band 5
		gc()
		Band_5_rad_vec <- getValues(Band_5_rad) # Vectors		
		gc()
		B5.OUT.PIX <- unique (c (B4.OUT.PIX, which (Band_5_rad_vec == Band_5_rad@data@min)))
		rm (Band_5_rad)
		rm (B4.OUT.PIX )
		rm (Band_5_rad_vec)
		gc()

		Band_6_rad <- raster (INPUT_BANDS [6]) # Band 6
		gc()
		Band_6_rad_vec <- getValues(Band_6_rad) # Vectors		
		gc()
		B6.OUT.PIX <- unique (c (B5.OUT.PIX, which (Band_6_rad_vec == Band_6_rad@data@min)))
		rm (Band_6_rad)
		rm (B5.OUT.PIX )
		rm (Band_6_rad_vec)
		gc()

		Band_7_rad <- raster (INPUT_BANDS [7]) # Band 7
		gc()
		Band_7_rad_vec <- getValues(Band_7_rad) # Vectors		
		gc()
		B7.OUT.PIX <- unique (c (B6.OUT.PIX, which (Band_7_rad_vec == Band_7_rad@data@min)))
		rm (Band_7_rad)
		rm (B6.OUT.PIX )
		rm (Band_7_rad_vec)
		gc()		
	
	# Digital numbers
		Band_1_DN <- raster (INPUT_BANDS_DN [1])[] # Band 1	
		B1Satu <- which (Band_1_DN == L7.FILLVALUE)  # Identify pixels that are saturated in the visible L7.PASS.CHOICE. Zhu has value of 255
		rm (Band_1_DN)
		gc()
		Band_2_DN <- raster (INPUT_BANDS_DN [2])[] # Band 2
		B2Satu <- which (Band_2_DN == L7.FILLVALUE)
		rm (Band_2_DN)
		gc()
		Band_3_DN <- raster (INPUT_BANDS_DN [3])[] # Band 3	
		B3Satu <- which (Band_3_DN == L7.FILLVALUE)	
		rm (Band_3_DN)
		gc()		
	## Clean
		B_Satu_unique_all <- sort (unique (c (B1Satu, B2Satu, B3Satu)))
		gc()
		B_Satu_table <- table( as.numeric (c (B1Satu, B2Satu, B3Satu)))
		dput (B1Satu, "B1Satu.RData"); dput (B2Satu, "B2Satu.RData"); dput (B3Satu, "B3Satu.RData");
		rm (B1Satu); rm(B2Satu); rm(B3Satu)
		gc()
		B_Satu_rows_all <- as.numeric (names (B_Satu_table[which (B_Satu_table == 3)]))
	 # Missing values
		id_missing <- B7.OUT.PIX 
		rm (B7.OUT.PIX)

	# Reflectance and temperature
		Band_1 <- raster (INPUT_BANDS [1]) # Band 1		
		gc()
		Band_1_vec_na <- Band_1[]
		rm (Band_1);  # Remove some of the rasters
		if (length (id_missing) > 0){ # Set NAs
			gc()
			Band_1_vec_na[id_missing] <- NA
		}
		dput(Band_1_vec_na, "Band_1_vec_na.RData")
		rm (Band_1_vec_na)

		Band_2 <- raster (INPUT_BANDS [2]) # Band 2		
		gc()
		Band_2_vec_na <- Band_2[]
		rm (Band_2);  # Remove some of the rasters
		if (length (id_missing) > 0){ # Set NAs
			gc()
			Band_2_vec_na[id_missing] <- NA
		}
		dput(Band_2_vec_na, "Band_2_vec_na.RData")
		rm (Band_2_vec_na)

		Band_3 <- raster (INPUT_BANDS [3]) # Band 3
		gc()
		Band_3_vec_na <- Band_3[]
		rm (Band_3);
		if (length (id_missing) > 0){ # Set NAs
			gc()
			Band_3_vec_na[id_missing] <- NA
		}
		dput(Band_3_vec_na, "Band_3_vec_na.RData")
		rm (Band_3_vec_na)
			  
		Band_4 <- raster (INPUT_BANDS [4]) # Band 4	
		gc()
		Band_4_vec_na <- Band_4[]
		rm (Band_4);	
		if (length (id_missing) > 0){ # Set NAs
			gc()
			Band_4_vec_na[id_missing] <- NA
		}	
		dput(Band_4_vec_na, "Band_4_vec_na.RData")
		rm (Band_4_vec_na)
			
		Band_5 <- raster (INPUT_BANDS [5]) # Band 5
		gc()
		Band_5_vec_na <- Band_5[]
		rm (Band_5);
		if (length (id_missing) > 0){ # Set NAs
			gc()
			Band_5_vec_na[id_missing] <- NA
		}	
		dput(Band_5_vec_na, "Band_5_vec_na.RData")
		rm (Band_5_vec_na)

		Band_6 <- raster (INPUT_BANDS [6]) # Band 6 in Celcius	
		gc()
		Band_6_vec_na <- Band_6[]
		rm (Band_6); 	
		if (length (id_missing) > 0){ # Set NAs
			gc()
			Band_6_vec_na[id_missing] <- NA
		}
		dput(Band_6_vec_na, "Band_6_vec_na.RData")
		rm (Band_6_vec_na)
			
		Band_7 <- raster (INPUT_BANDS [7]) # Band 7	
		gc()
		Band_7_vec_na <- Band_7[]
		rm (Band_7);	
		if (length (id_missing) > 0){ # Set NAs
			gc()
			Band_7_vec_na[id_missing] <- NA
		}
		dput(Band_7_vec_na, "Band_7_vec_na.RData")
		rm (Band_7_vec_na)
		
  ## Do		
	 # If all values are not NA
	  Band_7_vec_na <- dget( "Band_7_vec_na.RData")
	  if (length (id_missing) != length (Band_7_vec_na)){
		 rm (Band_7_vec_na)
		 Band_1_vec_na <- dget( "Band_1_vec_na.RData")
		 Band_2_vec_na <- dget( "Band_2_vec_na.RData")
		 Band_3_vec_na <- dget( "Band_3_vec_na.RData")
		 satu_Bv.data <- (Band_1_vec_na + Band_2_vec_na + Band_3_vec_na)[B_Satu_unique_all]
		 rm (Band_1_vec_na); rm (Band_2_vec_na); rm (Band_3_vec_na);
		 satu_Bv <- satu_Bv.data >= 1
		 rm ( satu_Bv.data) 
		 gc()
		 Satu_rows <- B_Satu_unique_all[which (satu_Bv == TRUE)]
		 rm (B_Satu_unique_all); rm (satu_Bv);
    # Potential cloud layer -- pass one
		Band_6_vec_na <- dget( "Band_6_vec_na.RData")
		BT <- Band_6_vec_na # Renamed (same as matlab name)
		rm (Band_6_vec_na);
			## Eq 1: Basic cloud test
				Band_3_vec_na <- dget( "Band_3_vec_na.RData");  Band_4_vec_na <- dget( "Band_4_vec_na.RData")
				NDVI <- (Band_4_vec_na - Band_3_vec_na) / (Band_4_vec_na + Band_3_vec_na) # Normalised Difference Vegetation Index
				dput(NDVI, "NDVI.RData")
				rm (Band_3_vec_na); rm (Band_4_vec_na); rm (NDVI);
				Band_2_vec_na <- dget( "Band_2_vec_na.RData");  Band_5_vec_na <- dget( "Band_5_vec_na.RData")
				NDSI <- (Band_2_vec_na - Band_5_vec_na) / (Band_2_vec_na + Band_5_vec_na) # Normalised Difference Snow Index	
				rm (Band_2_vec_na); rm (Band_5_vec_na);
				NDVI <- dget("NDVI.RData")		
				Band_7_vec_na <- dget( "Band_7_vec_na.RData")			
				Basic_Test <- ifelse (Band_7_vec_na > 300 & BT < 2700 & NDSI < 0.8 & NDVI < 0.8, TRUE, FALSE) # Run test
				dput (Basic_Test,"Basic_Test.RData");  rm (Basic_Test); rm (Band_7_vec_na); rm (BT)
				dput (NDSI,"NDSI.RData");  rm (NDSI);
				gc()
			## Eq 2: Whiteness test
				Band_1_vec_na <- dget( "Band_1_vec_na.RData"); Band_2_vec_na <- dget( "Band_2_vec_na.RData"); Band_3_vec_na <- dget( "Band_3_vec_na.RData");
				MeanVis <- (Band_1_vec_na + Band_2_vec_na + Band_3_vec_na) / 3
				B1.MeanVis.diff <- abs (Band_1_vec_na - MeanVis)
				rm (Band_1_vec_na)
				B2.MeanVis.diff <- abs (Band_2_vec_na - MeanVis)
				rm (Band_2_vec_na)
				B3.MeanVis.diff <- abs (Band_3_vec_na - MeanVis)
				rm (Band_3_vec_na)
				Bsum.MeanVis.diff <- (B1.MeanVis.diff + B2.MeanVis.diff + B3.MeanVis.diff) / MeanVis
				rm (MeanVis); rm (B1.MeanVis.diff); rm (B2.MeanVis.diff); rm (B3.MeanVis.diff);
				whiteness <- Bsum.MeanVis.diff
				rm (Bsum.MeanVis.diff)
				gc()
			## Adjust whiteness values
				# Set saturated whiteness to zero
				if (length (Satu_rows) > 0){
					whiteness[Satu_rows] <- 0
				}
				whiteness[whiteness > 1] <- 1 #C- whiteness goes over 1, but I don't know why. Therefore reset to 1. -C#	
				Whiteness_Test <- whiteness < 0.7
				rm(whiteness);
				dput(Whiteness_Test,"Whiteness_Test.RData")
				gc()
			## Eq 3: Haze test
				Band_1_vec_na <- dget( "Band_1_vec_na.RData"); Band_3_vec_na <- dget( "Band_3_vec_na.RData");
				HOT_Test.data <- Band_1_vec_na - 0.5 * Band_3_vec_na - 800 # Test values
				HOT_Test <- HOT_Test.data > 0 # Run test
				rm (Band_1_vec_na); rm (Band_3_vec_na); rm (HOT_Test.data)
				dput(HOT_Test,"HOT_Test.RData")
				rm (HOT_Test)
				gc()
			## Eq 4: Band 4 and 5 ratio
				Band_4_vec_na <- dget( "Band_4_vec_na.RData"); Band_5_vec_na <- dget( "Band_5_vec_na.RData");
				B45_Test <- (Band_4_vec_na / Band_5_vec_na) > 0.75 # Test
				rm (Band_4_vec_na); rm (Band_5_vec_na)
				dput (B45_Test, "B45_Test.RData")
				rm (B45_Test)
				gc()
			## Eq 5: Water test
				NDVI <- dget("NDVI.Rdata"); Band_4_vec_na <- dget( "Band_4_vec_na.RData");
				Water_Test <- ifelse ((NDVI < 0.01 & Band_4_vec_na < 1100)|(NDVI < 0.1 & NDVI > 0 & Band_4_vec_na < 500), TRUE, FALSE)
				rm (NDVI); rm (Band_4_vec_na);
				dput (Water_Test, "Water_Test.RData")
				rm (Water_Test)
				gc()
			## Eq 6: Potential cloud ?
				Basic_Test <- dget("Basic_Test.RData"); Whiteness_Test <- dget("Whiteness_Test.RData"); HOT_Test <- dget("HOT_Test.RData"); B45_Test <- dget("B45_Test.RData");
				PCP <- ifelse (Basic_Test == TRUE & Whiteness_Test == TRUE & HOT_Test == TRUE & B45_Test ==TRUE, TRUE, FALSE)
				rm (Basic_Test); rm (Whiteness_Test); rm (HOT_Test); rm (B45_Test);		
				gc()
			
			####S Choice switch
				choice_switch_1 <- length (which (PCP == TRUE)) / length (PCP)
				dput(PCP,"PCP.RData")
				rm (PCP)
				gc()
			####E
			
			if (choice_switch_1 >= 0.99) { # Too many clouds, go straight to end. Set entire tile to cloud.
				Band_1_vec_na <- dget( "Band_1_vec_na.RData"); 
				Band_1 <- raster (INPUT_BANDS [1]);
				Fmask <- rep (1, length (Band_1_vec_na)) # Set vector
				Fmask_ras <- setValues (Band_1, Fmask) # Turn vector into raster
				writeRaster(Fmask_ras,OUTPUT,overwrite = TRUE) # Save raster
				rm (Band_1_vec_na); rm (Band_1); rm (Fmask);rm (Fmask_ras);
			}	
			if (choice_switch_1 < 0.99){ # Not fully cloud. Potential cloud layer -- pass two 
				# Temperature probability for water
				if (run.prob.over.water == "yes"){ # Are there any bodies of water to run this over? 
					## Eq 7: Temperature test over water
						Water_Test <- dget("Water_Test.RData")
						Band_7_vec_na <- dget ("Band_7_vec_na.RData")
						Clearsky_Water <- (Water_Test == TRUE & Band_7_vec_na <= 300) #C- Inconsistency between matlab and paper? -C#
						rm (Water_Test); rm (Band_7_vec_na);
					## Eq 8
						Band_6_vec_na <- dget( "Band_6_vec_na.RData")
						BT <- Band_6_vec_na # Renamed (same as matlab name)
						rm (Band_6_vec_na);
						T_water <- quantile (as.numeric (BT[Clearsky_Water == TRUE]), 0.825, na.rm = TRUE)
						dput (T_water,"T_water.RData")
						gc()
					## Eq 9
						wTemperature_Prob <- (T_water - BT) / 400
						rm (BT)
						dput (wTemperature_Prob, "wTemperature_Prob.RData")
						rm (wTemperature_Prob)
						rm (T_water)
						gc ()
					## Eq 10: Brightness test over water
						t_bright <- 1100
						Band_5_vec_na <- dget ("Band_5_vec_na.RData")
						Brightness_Prob <- Band_5_vec_na / t_bright
						rm (Band_5_vec_na)
						Brightness_Prob[Brightness_Prob > 1] <- 1
					## Eq 11: Final probability mask over water
						T_water <- dget ("T_water.RData")
						if (is.na (T_water) == FALSE){
							rm (T_water)
							wTemperature_Prob <- dget ("wTemperature_Prob.RData")
							wCloud_Prob <- wTemperature_Prob * Brightness_Prob * 100
							rm (Brightness_Prob); rm (wTemperature_Prob)
							gc ()
							} else { #C- If no clear-sky water was identified -C#
							rm (T_water)
							wCloud_Prob <- 1 / (Brightness_Prob) 
							rm (Brightness_Prob)
							dput (wCloud_Prob,"wCloud_Prob.RData")
							gc()
						}
					}
			# Temperature probability for land
				## Eq 12
					T_buffer <- 4 * 100
					PCP <- dget ("PCP.RData"); Water_Test <- dget ("Water_Test.RData")
					Clearsky_Land <- ifelse (PCP == FALSE & Water_Test == FALSE, TRUE, FALSE)
					rm (PCP); rm (Water_Test);
					dput(Clearsky_Land,"Clearsky_Land.RData")
					gc()
				## Eq 13	
					Band_6_vec_na <- dget( "Band_6_vec_na.RData")
					BT <- Band_6_vec_na # Renamed (same as matlab name)
					rm (Band_6_vec_na);
					T_low <- quantile (as.numeric (na.omit (BT[Clearsky_Land == TRUE])), 0.175) # Low percentile
					T_high <- quantile (as.numeric (na.omit (BT[Clearsky_Land == TRUE])), 0.825) # High percentile
					rm (Clearsky_Land); rm (BT)
				# Other stuff
					t_tempL <- T_low - T_buffer
					t_tempH <- T_high + T_buffer
					rm (T_high); 
					dput(T_low, "T_low.RData"); rm (T_low);
					Temp_l <- t_tempH - t_tempL
					dput(t_tempL,"t_tempL.RData"); rm (t_tempL)
					
				## Eq 14
					Band_6_vec_na <- dget( "Band_6_vec_na.RData")
					BT <- Band_6_vec_na # Renamed (same as matlab name)
					rm (Band_6_vec_na);
					lTemperature_Prob <- (t_tempH - BT) / Temp_l
					rm (BT)
					dput(t_tempH,"t_tempH.RData"); rm (t_tempH)
					dput(Temp_l,"Temp_l.RData"); rm (Temp_l)
					lTemperature_Prob[lTemperature_Prob < 0] <- 0 #C- Probability can be > 1 -C#
				####S Choice switch
					Clearsky_Land <- dget("Clearsky_Land.RData")
					choice_switch_2 <- length (which (Clearsky_Land == TRUE)) / length (Clearsky_Land)
					rm (Clearsky_Land)
				####E
			}

				#C- Fmask will use the clear clear-sky pixels from both land and water for calculating temperature probability -C#
				if (choice_switch_1 < 0.99){
					if (choice_switch_2 < 0.001){ 	
					# NDSI
						NDSI <- dget ("NDSI.RData")
						modified_NDSI <- NDSI
						B3Satu <- dget ("B3Satu.RData")
						modified_NDSI_tf_1 <- NDSI[B3Satu]<0
						rm (NDSI)
						modified_NDSI_tf_2 <- B3Satu[which(modified_NDSI_tf_1==TRUE)]
						rm (B3Satu); rm (modified_NDSI_tf_1)
						modified_NDSI[modified_NDSI_tf_2] <- 0								
					## Eq 20: potential snow layer
						Band_4_vec_na <- dget ("Band_4_vec_na.RData"); Band_2_vec_na <- dget ("Band_2_vec_na.RData"); 
						Band_6_vec_na <- dget( "Band_6_vec_na.RData")
						BT <- Band_6_vec_na # Renamed (same as matlab name)
						rm (Band_6_vec_na);
						PSL <- ifelse(modified_NDSI>0.15&BT<380&Band_4_vec_na>1100&Band_2_vec_na>1000,TRUE,FALSE)
					# Vector with values 4) cloud pixels, 2) shadow pixels, 3) snow pixels
						Band_1 <- raster (INPUT_BANDS [1])
						Fmask <- rep(0,length(Band_1))
						rm (modified_NDSI); rm (Band_4_vec_na); rm (Band_2_vec_na); rm (BT);
						Fmask[PSL == TRUE] <- 4 	# Snow
						Fmask[id_missing] <- NA # NA post 2003
						Fmask_ras <- setValues(Band_1,Fmask) # Create raster
						writeRaster(Fmask_ras,OUTPUT,overwrite = TRUE) # Save raster
						rm (Band_1); rm (Fmask_ras);
					}
				}

			#C- Fmask will use the clear clear-sky pixels land only for calculating temperature probability -C#
				if (choice_switch_1 < 0.99){				
					if (choice_switch_2 >= 0.001){			
						## Variability probability
						 # Modify vales for saturated pixels
						  # NDVI
							NDVI <- dget ("NDVI.RData")
							modified_NDVI <- NDVI	
							B2Satu <- dget ("B2Satu.RData")
							modified_NDVI_tf_1 <- NDVI[B2Satu] > 0
							modified_NDVI_tf_2 <- B2Satu[which (modified_NDVI_tf_1 == TRUE)]
							rm(modified_NDVI_tf_1); rm (B2Satu) ; rm (NDVI)
							modified_NDVI[modified_NDVI_tf_2] <- 0
							dput(modified_NDVI,"modified_NDVI.RData")
							rm (modified_NDVI)
							# NDSI
								NDSI <- dget ("NDSI.RData")	
								modified_NDSI <- NDSI
								B3Satu <- dget ("B3Satu.RData")
								modified_NDSI_tf_1 <- NDSI[B3Satu] < 0
								rm (NDSI); 
								modified_NDSI_tf_2 <- B3Satu[which (modified_NDSI_tf_1==TRUE)]
								rm (B3Satu); rm (modified_NDSI_tf_1);
								modified_NDSI[modified_NDSI_tf_2] <- 0
								rm (modified_NDSI_tf_2);
								dput (modified_NDSI,"modified_NDSI.RData")
								rm (modified_NDSI);
							## Eq 15
#START###################################################
	# COMMENT
	# JH (25/4/12): Change from the Matlab or paper.					
							#C- I removed whiteness because it can equal 1. Hence, Variability prob would be 0. Therefore, subsequent probabilities would be 0 -C#
								modified_NDVI <- dget("modified_NDVI.RData"); modified_NDSI <- dget ("modified_NDSI.RData")
								Variability_Prob <- 1 #-max(abs(modified_NDVI),abs(modified_NDSI),na.rm=T)#,whiteness,na.rm=T)
								rm (modified_NDSI); rm (modified_NDVI)
								#Variability_Prob <- 1-max(abs(modified_NDVI))#,abs(modified_NDSI))
#END#####################################################							
						## Cloud probability for land
							## Eq 16
								lCloud_Prob <- (lTemperature_Prob * Variability_Prob) * 100
								rm (Variability_Prob)
							 ## Eq 17
								cldprob <- 22.5 # Default cloud probability threshold for cloud detection
								Clearsky_Land <- dget ("Clearsky_Land.RData")
								Land_threshold <- as.numeric (quantile (lCloud_Prob[Clearsky_Land == TRUE], 0.825, na.rm = TRUE) + cldprob) 
								rm (Clearsky_Land);
								# Treshold for defining cloud cover
					  ## Potential cloud layer
						 ## Eq 18
							PCP <- dget ("PCP.RData");  Water_Test <- dget ("Water_Test.RData");
							if (run.prob.over.water=="yes"){							
								clr_max <- 50 # Treshold
								wCloud_Prob <- dget ("wCloud_Prob.RData")
								PCL_t_1 <- PCP == TRUE & Water_Test == TRUE & wCloud_Prob > clr_max # Cloud over water
								rm (wCloud_Prob)
								PCL_t_1 <- dput (PCL_t_1, "PCL_t_1.RData")
								rm (PCL_t_1)
							}
 							PCL_t_2 <- PCP == TRUE & Water_Test == FALSE & lCloud_Prob > Land_threshold # Cloud over land
							rm (Land_threshold); rm (PCP);
							PCL_t_2.length <- length (PCL_t_2)
							PCL_t_2 <- dput (PCL_t_2, "PCL_t_2.RData")
							rm (PCL_t_2)
							PCL_t_3 <- lCloud_Prob > 99 & Water_Test == FALSE # High probability of cloud over land
							rm (Water_Test); rm (lCloud_Prob);
							PCL_t_3 <- dput (PCL_t_3, "PCL_t_3.RData")
							rm (PCL_t_3)
							BT <- dget( "Band_6_vec_na.RData") # Band6 and BT are the same thing, just change of nomenclature
							T_low <- dget ("T_low.RData")
							PCL_t_4 <- BT < (T_low - 3500) # Really cold cloud
							rm (BT); rm (T_low);
							PCL_t_4 <- dput (PCL_t_4, "PCL_t_4.RData")
							rm (PCL_t_4)
						## Classify each pixel
							PCL_raw_class <- rep (NA, PCL_t_2.length) # Blank vector
							rm (PCL_t_2.length)
							# Add values to vector
								if (run.prob.over.water=="yes"){
									PCL_t_1 <- dget ("PCL_t_1.RData")
									PCL_raw_class[PCL_t_1] <- 1 
									rm (PCL_t_1);
								}	
								PCL_t_2 <- dget ("PCL_t_2.RData")		
									PCL_raw_class[PCL_t_2] <- 2
								rm (PCL_t_2);
								PCL_t_3 <- dget ("PCL_t_3.RData")
									PCL_raw_class[PCL_t_3] <- 3
								rm (PCL_t_3);
								PCL_t_4 <- dget ("PCL_t_4.RData")
									PCL_raw_class[PCL_t_4] <- 4
								rm (PCL_t_4);
								PCL_raw_class[is.na(PCL_raw_class)] <- 0
							# Make a 0,1 mask		
								PCL <- PCL_raw_class
								rm (PCL_raw_class);
								PCL[PCL!=0] <- 1
								dput (PCL, "PCL.RData")
						## Cloud work
							# Clumping (group unique clouds)
								Band_1 <- raster (INPUT_BANDS [1])
								PCL_ras <- setValues(Band_1,PCL)
								writeRaster (PCL_ras, "PCL_ras.tif", overwrite = TRUE)
								rm (PCL); rm (Band_1);
								PCL_clump <- clump (PCL_ras)
								PCL_clump_vec <- PCL_clump[]
								PCL_clump.ncell <- ncell(PCL_clump)
								rm (PCL_clump)
							## Add 8 pixels around each cloud pixel in a clump
								PCL <- dget("PCL.RData")
								Cloud_pix_row <- which (PCL == 1) # Cloud pixels
								rm (PCL)
								dput (Cloud_pix_row, "Cloud_pix_row.RData")
						# Adjacency (outside the edge of cloud)
							PCL_ras <- raster ("PCL_ras.tif")
							PCL_adj_1 <- adjacency (PCL_ras, Cloud_pix_row, 1 : PCL_clump.ncell, directions = 8)
							rm (PCL_ras)
							rm (Cloud_pix_row);
							PCL_adj_v1 <- vector (length = PCL_clump.ncell)			
							PCL_adj_v1[PCL_adj_1[,2]] <- 2
							rm (PCL_adj_1)
							Cloud_pix_row <- dget ("Cloud_pix_row.RData")
							PCL_adj_v1[Cloud_pix_row] <- 1
							dput(PCL_adj_v1,"PCL_adj_v1.RData")
							rm (Cloud_pix_row)
						# Adjacency (at the edge of a cloud)
							Outer_edge_cloud_pix <- which (PCL_adj_v1 == 2)
							rm (PCL_adj_v1);
							PCL_ras <- raster ("PCL_ras.tif")
							PCL_adj_2 <- adjacency (PCL_ras, fromCells = Outer_edge_cloud_pix,toCells = 1:PCL_clump.ncell, directions = 8)
							rm (PCL_ras)
							PCL_adj_v2 <- vector(length = PCL_clump.ncell)	
							Cloud_pix_row <- dget ("Cloud_pix_row.RData")
							Inner_edge <- which (PCL_adj_2[,2] %in% Cloud_pix_row)
							rm (Cloud_pix_row)
							PCL_adj_v2[PCL_adj_2[,2]] <- 3
							rm (PCL_adj_2)
							PCL_adj_v2[PCL_adj_v2 != 3] <- 0
							PCL_adj_v2_2 <- rep (NA, length (PCL_adj_v2))
							PCL_adj_v1 <- dget ("PCL_adj_v1.RData")
							PCL_adj_v2_2[PCL_adj_v1==1&PCL_adj_v2==3] <- 1
							rm (PCL_adj_v2); rm (PCL_adj_v1)
							PCL_adj_v2_2[is.na(PCL_adj_v2_2)==TRUE] <- 0
							Cloud_edge_pix <- PCL_adj_v2_2
							rm (PCL_adj_v2_2)
						# Re-clump without edges
							PCL <- dget("PCL.RData")
							PCL_edge <- PCL
							rm (PCL)
							PCL_edge[which (Cloud_edge_pix == 1)] <- 0
							rm (Cloud_edge_pix)
							# Can't really clump if only 1 pixel
								Band_1 <- raster (INPUT_BANDS [1]) # Band 1		
								PCL_ras_edge <- setValues (Band_1,PCL_edge)
								if (length(PCL_edge[PCL_edge==1]) > 1){
									PCL_clump_edge <- clump (PCL_ras_edge)
									PCL_clump_vec_edge <- PCL_clump_edge[]
										} else {
									PCL_clump_edge <- PCL_ras_edge
									PCL_clump_vec_edge <- PCL_edge
								}
								rm (PCL_edge); rm (Band_1)
								dput(PCL_clump_vec_edge,"PCL_clump_vec_edge.RData")	
					## EQ 22/23 Cloud base 
						BT <- dget( "Band_6_vec_na.RData") # Band6 and BT are the same thing, just change of nomenclature
						BT_split <- split (BT, PCL_clump_vec_edge) # Split BT values
						rm (BT)
						# Length of clouds
							PCL_clump_vec_edge.table <- table (PCL_clump_vec_edge)
							PCL_clump_vec_edge_table_df <- data.frame (VAL = names (PCL_clump_vec_edge.table), COU = as.numeric (PCL_clump_vec_edge.table))
							rm (PCL_clump_vec_edge.table); rm (PCL_clump_vec_edge);
							Tcloud_base <- rep (NA, nrow (PCL_clump_vec_edge_table_df)) # Blank vector
					# Cloud base BT
						for(a in 1 : length(Tcloud_base)){
							print (paste("Cloud base",a,"of",length(Tcloud_base)))
							BT_split_iter <- BT_split[[a]]
							# Cloud_id <- PCL_clump_vec_edge_table_df[a,1] # Not used anywhere else. Don't know why here. (JH)
							R <- sqrt (PCL_clump_vec_edge_table_df[a,2] / (2 * pi))	
							if (R >= 8){ #C- Between 0 and 1 in R -C#
								quant.val <- (((R - 8)^2)) / R^2
								Tcloud_base[a] <- as.numeric (quantile (BT_split_iter,quant.val, na.rm = TRUE))
							}
							if (R < 8){
								Tcloud_base[a] <- min (BT_split_iter)
								Tcloud_base[a] <- ifelse (Tcloud_base[a] < 0, NA, Tcloud_base[a]) #c- read through more -c#
							}
						}
						dput(Tcloud_base,"Tcloud_base.RData")
						rm (PCL_clump_vec_edge_table_df); rm (BT_split_iter); rm (R);
					## EQ 23: Tcloud_object > base == base
						BT_split_mod <- BT_split
						for(a in 1:length (Tcloud_base)){
							print (paste("Cloud base EQ23",a,"of",length(Tcloud_base)))
							in_1 <- BT_split[[a]]
							in_2 <- Tcloud_base[a]
							in_1[in_1 > in_2] <- in_2
							BT_split_mod[[a]] <- in_1
							rm (in_1); rm (in_2);
						}
						rm (BT_split)
						PCL_clump_vec_edge <- dget ("PCL_clump_vec_edge.RData")
						BT_mod <- unsplit (BT_split_mod, PCL_clump_vec_edge) # Modified BT
						rm (BT_split_mod); rm (PCL_clump_vec_edge);
						dput(BT_mod,"BT_mod.RData")
					## Eq 20: potential snow layer
						modified_NDSI <- dget ("modified_NDSI.RData")
						BT <- dget( "Band_6_vec_na.RData") # Band6 and BT are the same thing, just change of nomenclature
						Band_4_vec_na <- dget( "Band_4_vec_na.RData")
						Band_2_vec_na <- dget( "Band_2_vec_na.RData")
						PSL <- ifelse (modified_NDSI > 0.15 & BT < 380 & Band_4_vec_na > 1100 & Band_2_vec_na > 1000, TRUE, FALSE)
						rm (modified_NDSI); rm (BT); rm (Band_4_vec_na); rm (Band_2_vec_na); 
						dput(PSL,"PSL.Rdata")
						rm (PSL)
					## Object-based cloud and cloud shadow match
						PCL <- dget("PCL.RData")
						revised_ptm <- 100 * (length (which (PCL == 1)) / length (PCL))
						rm (PCL);
						if (revised_ptm >= 90){ # Too many clouds.  Set entire layer to cloud.
							print("No shadow match because too much cloud")
							Band_1_vec_na <- dget( "Band_1_vec_na.RData")
							Band_1 <- raster (INPUT_BANDS [1])
							Fmask <- rep (1, length (Band_1_vec_na)) # Set vector values
							Fmask_ras <- setValues (Band_1, Fmask)
							rm (Band_1); rm (Band_1_vec_na);
							writeRaster (Fmask_ras, OUTPUT, overwrite = TRUE)
						}
						if (revised_ptm < 90){ # Not too many clouds
				# RSAGA
				# Remove spikes and dips
				# Band 4
					Band_1 <- raster (INPUT_BANDS [1])
					Band_4_vec_na <- dget("Band_4_vec_na.RData")
					Band_4_vec_ref <- setValues (Band_1, as.numeric (Band_4_vec_na))
					writeRaster(Band_4_vec_ref,"BAND_4_TEMP_filled.sgrd" , overwrite = TRUE)
					rm (Band_1); rm (Band_4_vec_na); rm (Band_4_vec_ref);
					rsaga.fill.sinks(paste(getwd(), "BAND_4_TEMP_filled.sgrd",sep="/"), paste(getwd(), "BAND_4_TEMP_filled.sgrd",sep="/"), method = "planchon.darboux.2001")
				# Band 5
					Band_1 <- raster (INPUT_BANDS [1])
					Band_5_vec_na <- dget("Band_5_vec_na.RData")
					Band_5_vec_ref <- setValues (Band_1, as.numeric (Band_5_vec_na))
					writeRaster(Band_5_vec_ref,  "BAND_5_TEMP.sgrd", overwrite = TRUE)
					rm (Band_1); rm (Band_5_vec_na); rm (Band_5_vec_ref);
					rsaga.fill.sinks(paste(getwd(), "BAND_5_TEMP_filled.sgrd",sep="/"), paste(getwd(), "BAND_5_TEMP_filled.sgrd",sep="/"), method = "planchon.darboux.2001")	
				# Re-format data	
					rsaga.sgrd.to.esri(in.sgrd=c("BAND_4_TEMP_filled.sgrd", "BAND_5_TEMP_filled.sgrd"), out.grids=c("BAND_4_TEMP_filled.tif", "BAND_5_TEMP_filled.tif"), out.path=getwd(), prec=1)			
				## Eq 19
					Band_4_filled_vec_na <- raster("BAND_4_TEMP_filled.tif")[]
					Band_5_filled_vec_na <- raster("BAND_5_TEMP_filled.tif")[]
					Band_4_vec_na <- dget("Band_4_vec_na.RData")
					Band_5_vec_na <- dget("Band_5_vec_na.RData")
					
					shadow_prob <- apply (cbind (Band_4_filled_vec_na - Band_4_vec_na,Band_5_filled_vec_na - Band_5_vec_na), 1, function(x) min (x, na.rm = TRUE))
					rm (Band_4_filled_vec_na); rm (Band_5_filled_vec_na); rm (Band_4_filled_vec); rm (Band_5_filled_vec);
#START###################################################
	# COMMENT
	# JH (28/9/12): Can't find PCSL used anywhere is fmask. Has to go somehwere?	
					PCSL <- shadow_prob > 200
#END#####################################################					
					rm (shadow_prob);
					rm (PCSL);
					## Constants
						t_obj <- Tcloud_base
						Tsimilar <- 0.3
						Tbuffer <- 0.98 # Treshold for matching buffering					
						num_cldoj <- 25 # Minimum matched cloud object (pixels)
						num_pix <- 8 # Number of inward pixels (240 M) for cloud base temperature
						rate_elapse <- 6.5 # Lapse rates
						rate_dlapse <- 9.8 # Lapse rates
						Max_cl_height <- 12000 	# Max cloud base height (m)
						Min_cl_height <- 200 # Min cloud base height (m)				
					## Cloud heights
						t_tempL <- dget ("t_tempL.RData")
						Min_cl_height_1 <- unlist (lapply (t_obj, function (x) max (Min_cl_height, 10 * (t_tempL - 400 - x) / rate_dlapse, na.rm = TRUE)))
						rm (t_tempL)
						Min_cl_height_1 <- ifelse (Min_cl_height_1 < 0, NA, Min_cl_height_1)
						t_tempH <- dget ("t_tempH.RData")
						Max_cl_height_1 <- unlist (lapply (t_obj, function (x) min (Max_cl_height, 10 * (t_tempH - 400 - x), na.rm = TRUE)))
						rm (t_tempH)
						Max_cl_height_1 <- ifelse (Max_cl_height_1 < 0, NA, Max_cl_height_1)
					## Object-based cloud and cloud shadow match
						BT_mod <- dget("BT_mod.RData")
						PCL_clump_vec_edge <- dget ("PCL_clump_vec_edge.RData")
						BT_mod_split <- split (BT_mod, PCL_clump_vec_edge) #C- Iteration for each cloud -C#
						rm (BT_mod)
						PCL_clump_vec_edge.mem <- data.frame (ID = 1:length (PCL_clump_vec_edge), PCL_clump_vec_edge = PCL_clump_vec_edge) #C- made due to memory problems -C#
						rm (PCL_clump_vec_edge)
						shadow_cells_list <- list()	
						h <- list()
						YX_mod_out_list <- list()
						## Each cloud
						Band_1 <-  raster (INPUT_BANDS [1])
						sub_size <- xres (Band_1) # Resolution (assumes xres==yres)
						rm (Band_1)
						
						
						
	
test.func.a <- function (){
	a = 1				
	h.raw <- 10 * (t_obj_a - temp_obj)	
	i_xy.raw <-  (sub_size * tan (L7.SUN.ELEVATION.TAN.RADIANS))
	X_move.raw 	<- cos (L7.SUN.AZIMUTH.TAN.RADIANS)
	Y_move.raw <-  sin (L7.SUN.AZIMUTH.TAN.RADIANS)
	b.list.1 <- vector("list", length(base_h))
	for(b in 1:2){#length (base_h)){
		#print (b)
		h <- h.raw / rate_elapse + base_h[b] # Each height in the cloud
		i_xy <- h / i_xy.raw 	# Shadow moved distance (pixel)
		# Movement (pixels)
		X_move <- round ((i_xy * X_move.raw) / 30)
		Y_move <- round ((i_xy * Y_move.raw) / 30)
		# Indentify cloud cells
		cloud.cell.id <- as.numeric (na.omit (PCL_clump_vec_edge.mem[PCL_clump_vec_edge.mem[,2] == a, 1]))
		YX_ori <- rowColFromCell (PCL_clump_edge, cloud.cell.id)	
		# Identify shadow rowCol positions
		YX_mod <- YX_ori
		out <- list(X_move,Y_move)
		#out <- list(YX_ori,YX_mod,X_move,Y_move)
		b.list.1 [[b]] <-  out
	}
	return (b.list.1)
	}
	system.time(test.func.a())
	Results.a <- test.func.a()
# user  system elapsed 
# 152.57    3.20  157.29 
	test.func.b <- function (){
	a = 1				
	h.raw <- 10 * (t_obj_a - temp_obj)	
	i_xy.raw <-  (sub_size * tan (L7.SUN.ELEVATION.TAN.RADIANS))
	X_move.raw 	<- cos (L7.SUN.AZIMUTH.TAN.RADIANS)
	Y_move.raw <-  sin (L7.SUN.AZIMUTH.TAN.RADIANS)
	b.list.1 <- vector("list", length(base_h))
	b.list.1 <- sapply(base_h[1:2], function (x) {
			#print(which(base_h == x))
			h <- h.raw / rate_elapse + x # Each height in the cloud
			i_xy <- h / i_xy.raw 	# Shadow moved distance (pixel)
			# Movement (pixels)
				X_move <- round ((i_xy * X_move.raw) / 30)
				Y_move <- round ((i_xy * Y_move.raw) / 30)
			# Indentify cloud cells
				cloud.cell.id <- as.numeric (na.omit (PCL_clump_vec_edge.mem[PCL_clump_vec_edge.mem[,2] == a, 1]))
				YX_ori <- rowColFromCell (PCL_clump_edge, cloud.cell.id)
			# Identify shadow rowCol positions
				YX_mod <- YX_ori
				#out <- list(YX_ori,YX_mod,X_move,Y_move)
				out <- list(X_move,Y_move)
			# Return
			return (out)
		})
	}
	system.time(test.func.b())
	Results.b <- test.func.b()
   # user  system elapsed 
# 1184.64    0.53 1187.20 

	test.func.c <- function (a){		
		h.raw <- 10 * (t_obj_a - temp_obj)	
		i_xy.raw <-  (sub_size * tan (L7.SUN.ELEVATION.TAN.RADIANS))
		X_move.raw 	<- cos (L7.SUN.AZIMUTH.TAN.RADIANS)
		Y_move.raw <-  sin (L7.SUN.AZIMUTH.TAN.RADIANS)
		#b.list.1 <- vector("list", length(base_h))
		b.list.1 <- sapply(base_h, function (x) {
				#print(paste(which(base_h == x), "of",length(base_h)))
				# Each height in the cloud
					h <- h.raw / rate_elapse + x
				return (h)
		})
		b.vec.1 <- b.list.1
		i_xy <- b.vec.1 / i_xy.raw
		# Movement (pixels)
			X_move <- (i_xy * X_move.raw) / 30
			Y_move <- (i_xy * Y_move.raw) / 30
		# Indentify cloud cells
			PCL_clump_vec_edge.mem.a <- PCL_clump_vec_edge.mem[,1]
			PCL_clump_vec_edge.mem.b <- PCL_clump_vec_edge.mem[,2]
			PCL_clump_vec_edge.mem.c <- which (PCL_clump_vec_edge.mem.b == a)
			PCL_clump_vec_edge.mem.d <- PCL_clump_vec_edge.mem.a[PCL_clump_vec_edge.mem.c]
			rm (PCL_clump_vec_edge.mem.a); rm (PCL_clump_vec_edge.mem.b);rm (PCL_clump_vec_edge.mem.c);
			cloud.cell.id <- as.numeric (na.omit (PCL_clump_vec_edge.mem.d))
			YX_ori <- rowColFromCell (PCL_clump_edge, cloud.cell.id)	
		# Identify shadow rowCol positions
			YX_mod <- vector("list", 2)
			length (YX_ori[, 1] - Y_move)
		if (L7.SUN.AZIMUTH < 180){
			YX_mod[[1]] <- YX_ori[,1] - Y_move
			YX_mod[[2]] <- YX_ori[,2] - X_move
			} else {
			YX_mod[[1]] <- YX_ori[,1] + Y_move
			YX_mod[[2]] <- YX_ori[,2] + X_move
			rm(YX_ori)
		}
		# Values that are off the map
			#YX_mod[YX_mod[,1] <= 0,] <- NA
			YX_mod[[1]][YX_mod[[1]][,] <= 0] <- NA
			YX_mod[[2]][YX_mod[[2]][,] <= 0] <- NA
			#YX_mod[[1]] <- apply(YX_mod[[1]],2,function (x) ifelse (x <0, NA, x))
			#YX_mod[[2]] <- apply(YX_mod[[2]],2,function (x) ifelse (x <0, NA, x))
		# Merge YX_mod
			o1 <- c(YX_mod[[1]][,1],YX_mod[[2]][,1])
			o2 <- matrix (o1,ncol=2)
			YX_mod_merge_1 <- na.omit(o2)
			attributes(YX_mod_merge_1)$na.action <- NULL
			rm (o1); rm (o2)
			YX_mod_merge_2 <- vector("list", (ncol(YX_mod[[1]])-1))
				for(i in 1:(ncol(YX_mod[[1]])-1)){
					o1 <- c(YX_mod[[1]][,(i+1)],YX_mod[[2]][,(i+1)])
					o2 <- matrix (o1,ncol=2)
					rm (o1)
					YX_mod_merge_2[[i]] <- na.omit(o2)
					attributes(YX_mod_merge_2[[i]])$na.action <- NULL
				}
		# Shadow cells	
			shadow_cells_list_out_b <- shad_cell_func (PCL_clump_edge, as.data.frame (YX_mod_merge_1, stringsAsFactors = FALSE))
			#count <- 0
			shadow_cells_list_out_b <-sapply(YX_mod_merge_2, function (x){
				#count <- ifelse(count!=0,count+1,1)
				#print (count)
				o1 <-  as.numeric (unique (na.omit (c (shadow_cells_list_out_b, shad_cell_func (PCL_clump_edge, as.data.frame (x, stringsAsFactors = FALSE))))))
				return (o1)
				})	
			return (shadow_cells_list_out_b)
		} 
		
#START###################################################
	# COMMENT
	# JH (28/9/12): Parallel computing may speed up this loop. It goes through each cloud, which means that iterations (loop of a) are independent.
		for(a in 1:length (Tcloud_base)){
			print(paste(a,"of",length (Tcloud_base)))
			if (is.na (Min_cl_height_1[a]) == FALSE & is.na (Max_cl_height_1[a]) == FALSE){
				shadow_cells_list_out_b <- test.func.c(a)
			} else {
				shadow_cells_list_out_b <- NA
			}
		}
	# JH (29/9/12): Mistmatch for some vector lenghts	
#END#####################################################
		dput(shadow_cells_list_out_b,"shadow_cells_list_out_b.RData")
	
						
						
						

						# for(a in 1:length (Tcloud_base)){
							# if (is.na (Min_cl_height_1[a]) == FALSE & is.na (Max_cl_height_1[a]) == FALSE){
								# shadow_cells_list_out <- list() # List
								# #i_step <- 2* # Move 2 pixels at a time
								# # Calculate cloud DEM with initial base height
	# #START###################################################
		# # COMMENT
		# # JH (28/9/12): Problem--minimum cloud height is bigger than max height)	
						# Min_cl_height_1[a] <- ifelse(Min_cl_height_1[a]>Max_cl_height_1[a],200,Min_cl_height_1) # Temporary work-around
						# base_h <- seq (Min_cl_height_1[a], Max_cl_height_1[a], 10)
	# #END#####################################################
								# t_obj_a <- t_obj[a]
								# #base_h <- Min_cl_height_1[a]
								# temp_obj <- BT_mod_split[[a]]
								# ## For each height of cloud
								# for(b in 1:length (base_h)){
									# h <- 10 * (t_obj_a-temp_obj) / rate_elapse + base_h[b] # Each height in the cloud
									# i_xy <- h / (sub_size * tan (L7.SUN.ELEVATION.TAN.RADIANS)) 	# Shadow moved distance (pixel)
									# # Movement (pixels)
									# X_move <- round ((i_xy * cos (L7.SUN.AZIMUTH.TAN.RADIANS)) / 30)
									# Y_move <- round ((i_xy * sin (L7.SUN.AZIMUTH.TAN.RADIANS)) / 30)
									# # Indentify cloud cells
									# cloud.cell.id <- as.numeric (na.omit (
									# PCL_clump_vec_edge.mem[PCL_clump_vec_edge.mem[,2] == a, 1]))
									# #which(PCL_clump_vec_edge==a)	
									# YX_ori <- rowColFromCell (PCL_clump_edge, cloud.cell.id)	
									# # Identify shadow rowCol positions
									# YX_mod <- YX_ori
								# if (L7.SUN.AZIMUTH < 180){
									# #if (L7.SUN.AZIMUTH.TAN.RADIANS < 180){
									# YX_mod[,1] <- YX_ori[,1] - Y_move
									# YX_mod[,2] <- YX_ori[,2] - X_move
									# } else {
									# YX_mod[,1] <- YX_ori[,1] + Y_move
									# YX_mod[,2] <- YX_ori[,2] + X_move
									# rm(YX_ori)
								# }
								# # Values that are off the map
								# YX_mod[YX_mod[,1] <= 0,] <- NA
								# YX_mod[YX_mod[,2] <= 0,] <- NA
								# YX_mod_out <- na.omit (YX_mod)
								# rm(YX_mod)
								# if (b == 1){
									# shadow_cells_list_out_b <- shad_cell_func (PCL_clump_edge, as.data.frame (YX_mod_out, stringsAsFactors = FALSE))
									# } else {
									# shadow_cells_list_out_b <- as.numeric (unique (na.omit (c (shadow_cells_list_out_b, shad_cell_func (PCL_clump_edge, as.data.frame (YX_mod_out, stringsAsFactors = FALSE))))))
									# rm(YX_mod_out)
								# }												
								# print(paste(b, length(base_h),sep=" of ")) # Print progress
								# }
								# } else {
									# shadow_cells_list_out_b <- NA
								# }
								# # Output for each cloud				
								# if (a == 1){
									# shadow_cells_list <- shadow_cells_list_out_b
									# rm(shadow_cells_list_out_b)
								# } else {
									# shadow_cells_list <- as.numeric (unique (na.omit (c (shadow_cells_list, shadow_cells_list_out_b))))
									# rm(shadow_cells_list_out_b)
								# }
								# # See where it is at
									# print("#####################################")
									# print (paste (a, length (Tcloud_base), sep = " of "))
					 # }
# #END#####################################################					 
					## Cloud raster
					Unique.cloud.pix <- unique(shadow_cells_list) # Unique pixels only
					if (length(Unique.cloud.pix)!=0){
						cloud_shadow <- rep (NA, length (PCL_clump_edge[])) # Make vector
						cloud_shadow[unlist (Unique.cloud.pix)] <- 1 # Shadow	
						Band_1 <-  raster (INPUT_BANDS [1])						
						cloud_shadow_ras_1 <- setValues (Band_1, cloud_shadow) # Make raster
						rm (cloud_shadow)
						## Dialte cloud and shadow pixels 3 times. Zhu and Woodcock, 2012 call for buffers by 3 pixels in 8-connected directions -C#
						# 1 pixel out
							shadow_cal_1 <- adjacency (cloud_shadow_ras_1, fromCells = Unique.cloud.pix, toCells = c(1:ncell (cloud_shadow_ras_1)), directions = 8)
							rm (cloud_shadow_ras_1)
							shadow_cal_1_vec <- rep (NA, length (PCL_clump_edge[]))		
							shadow_cal_1_vec[shadow_cal_1[,2]] <- 1
							shadow_cal_1_ras <- setValues (Band_1, shadow_cal_1_vec)				
						# 2 pixels out
							shadow_cal_2_vec_vals <- which (shadow_cal_1_vec == 1)
							rm (shadow_cal_1_vec)
							shadow_cal_2 <- adjacency (shadow_cal_1_ras, fromCells = shadow_cal_2_vec_vals, toCells = c(1:ncell (shadow_cal_1)), directions = 8)
							rm (shadow_cal_1_ras)
							shadow_cal_2_vec <- rep (NA, length (PCL_clump_edge[]))		
							shadow_cal_2_vec[shadow_cal_2[,2]] <- 1
							rm (shadow_cal_2)
							shadow_cal_2_ras <- setValues (Band_1, shadow_cal_2_vec)
						# 3 pixels out
							shadow_cal_3_vec_vals <- which (shadow_cal_2_vec == 1)
							rm (shadow_cal_2_vec)
							shadow_cal_3 <- adjacency (shadow_cal_2_ras, fromCells = shadow_cal_3_vec_vals, toCells = c(1:ncell (shadow_cal_1)), directions = 8)
							rm (shadow_cal_2_ras)
							rm(shadow_cal_1)
							shadow_cal_3_vec <- rep (NA, length (PCL_clump_edge[]))		
							shadow_cal_3_vec[shadow_cal_3[,2]] <- 1
						} else {
							shadow_cal_3_vec <- rep (0, length (PCL_clump_edge[]))
						}
						
					## Classify output #C- Priority: 1) cloud pixels, 2) shadow pixels, 3) snow pixels
						Fmask <- rep(0,length(PCL_clump_edge[]))						
						PSL <- dget("PSL.RData")
						Fmask[PSL==TRUE] <- 3 # Snow
						rm (PSL)
						Fmask[shadow_cal_3_vec[]==1] <- 2 # Shadows
						rm (shadow_cal_3_vec)
						PCL <- dget("PCL.RData")
						Fmask[PCL==1] <- 1 # Cloud
						rm (PCL)
						Fmask[id_missing] <- NA # NA post 2003
						Fmask_ras <- setValues(Band_1,Fmask) # Raster
						writeRaster(Fmask_ras, OUTPUT, overwrite = TRUE)
					}
				}
			 }
	    }
}

